#include "FDK.h"
#include "../Common/ConfigFile.h"
#include "../Common/CommonDef.h"

#include "math.h"
#include "time.h"
#include "assert.h"

#include <iostream>
#include <fstream>
#include <string>
using namespace std;

#include "ipp_w7.h"
#include "ipp.h"

#pragma comment( lib, "ippsmerged" )
#pragma comment( lib, "ippmmerged" )

FDKPara para;
float * pfBuffer = NULL;
float * pfConvKernal = NULL;
float * pfProjection = NULL;
float * pfImage = NULL;
float * pfParkerWeighting = NULL;
float * pfPmartices = NULL;
bool bPmatricesUsed = false;

int main(int argc, char *argv[])
{
	CConfigFile configfile;

	if(argc != 2 && argc != 4 && argc != 5)
	{
		cout << "Usage:\n FDK ConfigFilename [ProjectionFilename ReconstructionFilename [PmatricesFilename]]\n";
		return -1;
	}

	if(configfile.Open(argv[1]) == -1)
	{
		cout << "Failed to open the config file.\n";
		return -1;
	}

	char buf[256];
	if(argc == 2)
	{
		PARA_GET_STRING(configfile,"File","ProjectionFilename",para.strProjectionFilename,buf);
		PARA_GET_STRING(configfile,"File","ReconstructionFilename",para.strReconstructionFilename,buf);
		PARA_GET_STRING(configfile,"File","PmatricesFilename",para.strPmatricesFilename,buf);
	}
	else if(argc >= 4)
	{
		para.strProjectionFilename = argv[2];
		para.strReconstructionFilename = argv[3];
		if(argc == 5)
			para.strPmatricesFilename = argv[4];
	}
	
	PARA_GET_INT(configfile,"Scan","SourceType",para.nSourceType,buf);
	PARA_GET_FLOAT(configfile,"Scan","SourceDetectorDistance",para.fSourceDetectorDistance,buf);
	PARA_GET_FLOAT(configfile,"Scan","SourceAxisDistance",para.fSourceAxisDistance,buf);
	PARA_GET_FLOAT(configfile,"Scan","StartOffAngle",para.fStartOffAngle,buf);
	PARA_GET_FLOAT(configfile,"Scan","StartSourceHeight",para.fStartSourceHeight,buf);
	PARA_GET_INT(configfile,"Scan","ProjectionsPerRotation",para.nProjectionsPerRotation,buf);
	PARA_GET_INT(configfile,"Scan","ProjectionNumber",para.nProjectionNumber,buf);
	PARA_GET_INT(configfile,"Scan","DetectorType",para.nDetectorKind,buf);
	PARA_GET_FLOAT(configfile,"Scan","DetectorTotalWidth",para.fDetectorTotalWidth,buf);
	PARA_GET_FLOAT(configfile,"Scan","DetectorTotalHeight",para.fDetectorTotalHeight,buf);
	PARA_GET_INT(configfile,"Scan","DetectorColumnNumber",para.nDetectorColumnNumber,buf);
	PARA_GET_INT(configfile,"Scan","DetectorRowNumber",para.nDetectorRowNumber,buf);
	PARA_GET_FLOAT(configfile,"Scan","DetectorOffset",para.fOffSet,buf);
	PARA_GET_INT(configfile,"Image","PixelNumberX",para.nImageX,buf);
	PARA_GET_INT(configfile,"Image","PixelNumberY",para.nImageY,buf);
	PARA_GET_INT(configfile,"Image","PixelNumberZ",para.nImageZ,buf);
	PARA_GET_FLOAT(configfile,"Image","CenterX",para.fCenterX,buf);
	PARA_GET_FLOAT(configfile,"Image","CenterY",para.fCenterY,buf);
	PARA_GET_FLOAT(configfile,"Image","CenterZ",para.fCenterZ,buf);
	PARA_GET_FLOAT(configfile,"Image","SampleIntervalX",para.fSampleIntervalX,buf);
	PARA_GET_FLOAT(configfile,"Image","SampleIntervalY",para.fSampleIntervalY,buf);
	PARA_GET_FLOAT(configfile,"Image","SampleIntervalZ",para.fSampleIntervalZ,buf);
	PARA_GET_INT(configfile,"Reconstruction","FilterType",para.nFilterType,buf);
	PARA_GET_BOOL(configfile,"Reconstruction","ParkWeighting",para.bParkWeighting,buf);
	PARA_GET_BOOL(configfile,"Reconstruction","ReconstructCylinder",para.bReconstructCylinder,buf);
	
	InitializeGlobalData();
	FDK();
	CleanupGlobalData();

	return 0;
}

void InitializeGlobalData()
{
	float fDetectorWidthReciprocal;
	int nIndexRow;

	pfConvKernal = new float[2*para.nDetectorColumnNumber-1];
 	memset(pfConvKernal,0,(2*para.nDetectorColumnNumber-1)*sizeof(float));

	if(para.nDetectorKind == 0)
	{
		fDetectorWidthReciprocal = para.nDetectorColumnNumber/(para.fDetectorTotalWidth*PI/180.0f);
		if(para.nFilterType==1)
		{
			//////////////////////////////////////////////////////////////////////////
			//  S-L kernal
		for(nIndexRow=-para.nDetectorColumnNumber+1;nIndexRow<=para.nDetectorColumnNumber-1;nIndexRow++)
			pfConvKernal[nIndexRow+para.nDetectorColumnNumber-1]= float(
					-(para.fDetectorTotalWidth/(para.nDetectorColumnNumber*sin(para.fDetectorTotalWidth/para.nDetectorColumnNumber)))
					*(para.fDetectorTotalWidth/(para.nDetectorColumnNumber*sin(para.fDetectorTotalWidth/para.nDetectorColumnNumber)))
					*2*fDetectorWidthReciprocal/(para.fSourceAxisDistance*(4*nIndexRow*nIndexRow-1)*PI*PI));
		}
		else if(para.nFilterType==0)
		{
			//////////////////////////////////////////////////////////////////////////
			//  ramp kernal
			pfConvKernal[para.nDetectorColumnNumber-1]=float(0.25*fDetectorWidthReciprocal);
			for(nIndexRow=-para.nDetectorColumnNumber+1;nIndexRow<=para.nDetectorColumnNumber-1;nIndexRow++)
			{
				if (nIndexRow%2!=0)
					pfConvKernal[nIndexRow+para.nDetectorColumnNumber-1]=-fDetectorWidthReciprocal/(nIndexRow*nIndexRow*PI*PI);
			}
		}
		else
		{
			//////////////////////////////////////////////////////////////////////////
			//  Delta conv kernel
			pfConvKernal[para.nDetectorColumnNumber-1]=1;
		}
	}
	else if(para.nDetectorKind == 1)
	{
		fDetectorWidthReciprocal = para.nDetectorColumnNumber/para.fDetectorTotalWidth;
		if(para.nFilterType==1)
		{
			//////////////////////////////////////////////////////////////////////////
			//  S-L kernal
			for(nIndexRow=-para.nDetectorColumnNumber+1;nIndexRow<=para.nDetectorColumnNumber-1;nIndexRow++)
				pfConvKernal[nIndexRow+para.nDetectorColumnNumber-1]=-2*fDetectorWidthReciprocal/((4*nIndexRow*nIndexRow-1)*PI*PI);
		}
		else if(para.nFilterType==0)
		{
			//////////////////////////////////////////////////////////////////////////
			//  ramp kernal
			pfConvKernal[para.nDetectorColumnNumber-1]=float(0.25*fDetectorWidthReciprocal);
			for(nIndexRow=-para.nDetectorColumnNumber+1;nIndexRow<=para.nDetectorColumnNumber-1;nIndexRow++)
			{
				if (nIndexRow%2!=0)
					pfConvKernal[nIndexRow+para.nDetectorColumnNumber-1]=-fDetectorWidthReciprocal/(nIndexRow*nIndexRow*PI*PI);
			}
		}
		else
		{
			//////////////////////////////////////////////////////////////////////////
			//  Delta conv kernel
			pfConvKernal[para.nDetectorColumnNumber-1]=1;
		}
	}
	else
	{
		cout << "Unkonw Detector Type:" << para.nDetectorKind << endl;
	}

	//  Modify the conv kernal for equi-angular detector
	float fDetectorX;
	if(para.nDetectorKind == 0)
	{
		for(nIndexRow=-para.nDetectorColumnNumber+1;nIndexRow<=para.nDetectorColumnNumber-1;nIndexRow++)
		{
			fDetectorX=nIndexRow*para.fDetectorTotalWidth*PI/(para.nDetectorColumnNumber*180);
			if (nIndexRow)
				pfConvKernal[nIndexRow+para.nDetectorColumnNumber-1]*=fDetectorX*fDetectorX/(sin(fDetectorX)*sin(fDetectorX));
		}
	}
	
	pfBuffer = new float[3*para.nDetectorColumnNumber-2];
	pfImage = new float[para.nImageX*para.nImageY*para.nImageZ];
	memset(pfImage,0,para.nImageX*para.nImageY*para.nImageZ*sizeof(float));
	pfProjection = new float[para.nDetectorColumnNumber*para.nDetectorRowNumber];

	if(para.bParkWeighting)
	{
		cout << "Parker Weighting applied" << endl;
		pfParkerWeighting = new float[para.nDetectorColumnNumber*para.nProjectionNumber];
		memset(pfParkerWeighting, 0, para.nDetectorColumnNumber*para.nProjectionNumber*sizeof(float));
		int na,nc;
		float * p = pfParkerWeighting;
		float Ym,Y,t;
		float t0,t1,t2,t3,tmp;

		if(para.nDetectorKind == 0)	// Cylinder detector
			Ym = para.fDetectorTotalWidth*PI/360;
		else						// Flat Panel detector
			Ym = (float)atan(para.fDetectorTotalWidth*0.5f/para.fSourceDetectorDistance);

		Ym *= float(para.nProjectionNumber-1)/float(para.nProjectionNumber);
		t0 = -PI/2-Ym;
		t3 = PI/2+Ym;

		for(na=0;na<para.nProjectionNumber;na++)
		{
			t = (na-(para.nProjectionNumber-1)*0.5f)*2*PI/para.nProjectionsPerRotation;
			for(nc=0;nc<para.nDetectorColumnNumber;nc++)
			{
				tmp = (nc+0.5f)*para.fDetectorTotalWidth/para.nDetectorColumnNumber-0.5f*para.fDetectorTotalWidth;

				if(para.nDetectorKind == 0)	// Cylinder detector
					Y = tmp*PI/180;
				else						// Flat Panel detector
					Y = (float)atan(tmp/para.fSourceDetectorDistance);
				
				t1 = -PI/2+Ym-2*Y;
				t2 = PI/2-Ym-2*Y;

				if(t0 <= t && t < t1)
				{
					tmp = (float)sin(PI*0.25f*(t+PI/2+Ym)/(Ym-Y));
					*p = tmp*tmp*2;
//					tmp = 0.5f*(t+PI/2+Ym)/(Ym-Y);
//					*p = tmp*tmp*(3-2*tmp)*2;
				}
				else if (t1 <= t && t< t2)
				{
					*p = 2.0f;
				}
				else if (t2 <= t && t< t3)
				{
					tmp = (float)sin(PI*0.25f*(-t+PI/2+Ym)/(Ym+Y));
					*p = tmp*tmp*2;
//					tmp = 0.5f*(-t+PI/2+Ym)/(Ym+Y);
//					*p = tmp*tmp*(3-2*tmp)*2;
				}
				else
				{
					*p = 0;
				}
				p++;
			}
		}
	}


	if(para.bReconstructCylinder)
	{
		cout << "Reconstruct cylinder FOV only." << endl;
	}

	if(!para.strPmatricesFilename.empty())
	{
		double * buf = new double[12*para.nProjectionNumber];
		int nlength,n1,n2,n3;
		ifstream in;
		in.open(para.strPmatricesFilename.c_str(),ios::in|ios::binary);
		if(in)
		{
			in.seekg(0, ios::end);
			nlength = in.tellg();
			if(nlength == 12*para.nProjectionNumber*sizeof(double))
			{
				in.seekg(0, ios::beg);
				in.read((char*)buf,nlength);
				pfPmartices = new float[12*para.nProjectionNumber];
				for(n1=0; n1<para.nProjectionNumber; n1++)
					for(n2=0; n2<4; n2++)
						for(n3=0; n3<3; n3++)
							pfPmartices[n1*12+n3*4+n2]=buf[n1*12+n2*3+n3];
				cout << "PMatrices file used." << endl;
				bPmatricesUsed = true;
			}
			else
			{
				cout << "The length of PMatrices file (" << nlength << ") is not consist with projection number (" 
					<< para.nProjectionNumber << ")." << endl;
			}
			in.close();
		}
		else
		{
			cout << "Failed to open PMatrices File:" << para.strPmatricesFilename << endl;
		}
		SAFE_DELETE(buf);
	}
}

void CleanupGlobalData()
{
	SAFE_DELETE(pfBuffer);
	SAFE_DELETE(pfConvKernal);
	SAFE_DELETE(pfProjection);
	SAFE_DELETE(pfImage);
	SAFE_DELETE(pfParkerWeighting);
}

void FDK()
{
	time_t ltime;
	time( &ltime );
	cout << "Enter FDK reconstruction at " << ctime(&ltime);

	string str;
	str = para.strProjectionFilename.substr(para.strProjectionFilename.length()-4,para.strProjectionFilename.length()-1);
	ifstream in1, in2;
	ofstream out;
	float * pfAnglePos = new float[para.nProjectionNumber];
	int nAngle;

	if(str.compare(".prj") == 0)
	{
		in1.open(para.strProjectionFilename.c_str(), ios::in|ios::binary);
		if(!in1)
		{
			cout << "Failed to open projection file: " << para.strProjectionFilename << endl;
			exit(-1);
		}

		string str = para.strProjectionFilename + ".AnglePos";

		in2.open(str.c_str(),ios::in | ios::binary);
		if(in2)
		{
			in2.read((char*)pfAnglePos, para.nProjectionNumber*sizeof(float));
			in2.close();

			cout << "Angle position file found.\n";

			for(nAngle=0;nAngle<para.nProjectionNumber;nAngle++)
			{
				pfAnglePos[nAngle] = pfAnglePos[nAngle]*PI/180+PI;
			}
		}
		else
		{
			for(nAngle=0;nAngle<para.nProjectionNumber;nAngle++)
			{
				pfAnglePos[nAngle] = para.fStartOffAngle*PI/180+nAngle*2*PI/para.nProjectionsPerRotation;
			}
		}

		for(nAngle=0;nAngle<para.nProjectionNumber;nAngle++)
		{
			cout << "Processing " << nAngle+1 << "/" <<para.nProjectionNumber << " angle\r";
			in1.read((char*)pfProjection,para.nDetectorColumnNumber*para.nDetectorRowNumber*sizeof(float));
			
			if(para.bParkWeighting)
				ApplyParkerWeighting(pfProjection,nAngle);
			
			if(bPmatricesUsed)
			{
				if(para.nDetectorKind == 0)
					ComputeOneAngleEA(pfProjection,pfImage,pfAnglePos[nAngle]);
				else
					ComputeOneAngleEDPmatrices(pfProjection,pfImage,pfPmartices+12*nAngle,pfAnglePos[nAngle]);
			}
			else
			{
				if(para.nDetectorKind == 0)
					ComputeOneAngleEA(pfProjection,pfImage,pfAnglePos[nAngle]);
				else
					ComputeOneAngleED(pfProjection,pfImage,pfAnglePos[nAngle]);
			}
		}
		time( &ltime );
		cout << "Reconstruction finished at " << ctime(&ltime);
		in1.close();
	}
	else if(str.compare(".key") == 0)
	{
	}

	cout << "Writing results... ";
	out.open(para.strReconstructionFilename.c_str(),ios::out|ios::binary);
	out.write((char*)pfImage,para.nImageX*para.nImageY*para.nImageZ*sizeof(float));
	out.close();
	cout << "Done." << endl;

	SAFE_DELETE(pfAnglePos);
}

void ApplyParkerWeighting(float *pPrj, int nAngle)
{
	for(int nRow=0; nRow<para.nDetectorRowNumber; nRow++)
		ippsMul_32f_I(pfParkerWeighting+nAngle*para.nDetectorColumnNumber,pPrj+nRow*para.nDetectorColumnNumber,para.nDetectorColumnNumber);
}


void ComputeOneAngleED(float *pPrj, float *pRcn, float fAngle)
{
	//////////////////////////////////////////////////////////////////////////
	//  filtering
	int nIndex=0;
	int nRow,nCol;
	float fDetectorX,fDetectorY;
	
	for(nRow=0;nRow<para.nDetectorRowNumber;nRow++)
	{
		nIndex=nRow*para.nDetectorColumnNumber;
		fDetectorY=(nRow+0.5f)*para.fDetectorTotalHeight/para.nDetectorRowNumber-0.5f*para.fDetectorTotalHeight;
		//////////////////////////////////////////////////////////////////////////
		//  Preweighting projection data
		for(nCol=0;nCol<para.nDetectorColumnNumber;nCol++)
		{
			fDetectorX=(nCol+0.5f)*para.fDetectorTotalWidth/para.nDetectorColumnNumber-0.5f*para.fDetectorTotalWidth;
			pPrj[nIndex+nCol]=pPrj[nIndex+nCol]*para.fSourceDetectorDistance
				/float(sqrt(fDetectorX*fDetectorX+fDetectorY*fDetectorY+para.fSourceDetectorDistance*para.fSourceDetectorDistance));
		}
		//////////////////////////////////////////////////////////////////////////
		//  convolution
		ippsConv_32f(pPrj+nIndex,para.nDetectorColumnNumber,pfConvKernal,para.nDetectorColumnNumber*2-1,pfBuffer);	
		memcpy(pPrj+nIndex,pfBuffer+para.nDetectorColumnNumber-1,para.nDetectorColumnNumber*sizeof(float));
	}
	
	//////////////////////////////////////////////////////////////////////////
	//  backprojection
	const float fdu=para.fDetectorTotalWidth/para.nDetectorColumnNumber;
	const float fdv=para.fDetectorTotalHeight/para.nDetectorRowNumber;
	const float fD=para.fSourceDetectorDistance;
	const float fR=para.fSourceAxisDistance;

	float fX,fY,fZ;	
	int nx,ny,nz;

	float fnu,fnv;
	float fSinAngle,fCosAngle;
	float fdeltaz;

	float fdeltau1,fdeltau2,fdeltav1,fdeltav2;
	float fui1,fui2,fuvi;
	float fr,fU,fWeight;

	const float fCoefficient=fD*fR*PI/para.nProjectionsPerRotation;
	const float fFovRadius=min(para.nImageX*para.fSampleIntervalX*0.5f,para.nImageY*para.fSampleIntervalY*0.5f);
	
	int nuindex,nvindex;
	nIndex=0;

	fSinAngle=(float)sin(fAngle);
	fCosAngle=(float)cos(fAngle);
	for(ny=0;ny<para.nImageY;ny++)
	{
		fY=para.fCenterY+(ny-0.5f*para.nImageY+0.5f)*para.fSampleIntervalY;
		for(nx=0;nx<para.nImageX;nx++)
		{
			fX=para.fCenterX+(nx-0.5f*para.nImageX+0.5f)*para.fSampleIntervalX;
			fr=float(sqrt(fX*fX+fY*fY));

			if(para.bReconstructCylinder && fr > fFovRadius)
				continue;
			
			fU=fR+fX*fCosAngle+fY*fSinAngle;
			fnu=fD*(-fX*fSinAngle+fY*fCosAngle)/(fU*fdu)+para.fOffSet/fdu;
			
			nuindex=int(fnu+0.5f*(para.nDetectorColumnNumber-1));
			if(nuindex>0 && nuindex<para.nDetectorColumnNumber-1)
			{
				fdeltau1=fnu+(para.nDetectorColumnNumber-1)*0.5f-nuindex;
				fdeltau2=1.0f-fdeltau1;
				fWeight=fCoefficient/(fU*fU);
				for(nz=0;nz<para.nImageZ;nz++)
				{
					fZ=para.fCenterZ+(nz-0.5f*para.nImageZ+0.5f)*para.fSampleIntervalZ;
					fdeltaz=fZ-para.fStartSourceHeight;
					fnv=fD*fdeltaz/(fU*fdv);
					nvindex=int(fnv+0.5f*para.nDetectorRowNumber);
					if(nvindex>=0 && nvindex<para.nDetectorRowNumber-1)
					{
						fdeltav1=fnv+para.nDetectorRowNumber-int(fnv+para.nDetectorRowNumber);
						fdeltav2=1.0f-fdeltav1;
						fui1=pfProjection[nvindex*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[nvindex*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						fui2=pfProjection[(nvindex+1)*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[(nvindex+1)*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						
						fuvi=fui1*fdeltav2+fui2*fdeltav1;
						pfImage[nz*para.nImageX*para.nImageY+ny*para.nImageX+nx]+=fWeight*fuvi;
					}
					else if(nvindex==para.nDetectorRowNumber-1)
					{
						fui1=pfProjection[nvindex*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[nvindex*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						pfImage[nz*para.nImageX*para.nImageY+ny*para.nImageX+nx]+=fWeight*fui1;
					}
				}
			}

		}
	}
}

void ComputeOneAngleEDPmatrices(float *pPrj, float *pRcn, float * pfPmartix, float fAngle)
{
	//////////////////////////////////////////////////////////////////////////
	//  filtering
	int nIndex=0;
	int nRow,nCol;
	float fDetectorX,fDetectorY;
	
	for(nRow=0;nRow<para.nDetectorRowNumber;nRow++)
	{
		nIndex=nRow*para.nDetectorColumnNumber;
		fDetectorY=(nRow+0.5f)*para.fDetectorTotalHeight/para.nDetectorRowNumber-0.5f*para.fDetectorTotalHeight;
		//////////////////////////////////////////////////////////////////////////
		//  Preweighting projection data
		for(nCol=0;nCol<para.nDetectorColumnNumber;nCol++)
		{
			fDetectorX=(nCol+0.5f)*para.fDetectorTotalWidth/para.nDetectorColumnNumber-0.5f*para.fDetectorTotalWidth;
			pPrj[nIndex+nCol]=pPrj[nIndex+nCol]*para.fSourceDetectorDistance
				/float(sqrt(fDetectorX*fDetectorX+fDetectorY*fDetectorY+para.fSourceDetectorDistance*para.fSourceDetectorDistance));
		}
		//////////////////////////////////////////////////////////////////////////
		//  convolution
		ippsConv_32f(pPrj+nIndex,para.nDetectorColumnNumber,pfConvKernal,para.nDetectorColumnNumber*2-1,pfBuffer);	
		memcpy(pPrj+nIndex,pfBuffer+para.nDetectorColumnNumber-1,para.nDetectorColumnNumber*sizeof(float));
	}
	
	//////////////////////////////////////////////////////////////////////////
	//  backprojection
	const float fdu=para.fDetectorTotalWidth/para.nDetectorColumnNumber;
	const float fdv=para.fDetectorTotalHeight/para.nDetectorRowNumber;
	const float fD=para.fSourceDetectorDistance;
	const float fR=para.fSourceAxisDistance;

	float fX,fY,fZ;	
	float fPos[4];
	float fResult[3];

	int nx,ny,nz;

	float fnu,fnv;
	float fSinAngle,fCosAngle;

	float fdeltau1,fdeltau2,fdeltav1,fdeltav2;
	float fui1,fui2,fuvi;
	float fr,fU,fWeight;
	

	const float fCoefficient=fD*fR*PI/para.nProjectionsPerRotation;
	const float fFovRadius=min(para.nImageX*para.fSampleIntervalX*0.5f,para.nImageY*para.fSampleIntervalY*0.5f);
	
	int nuindex,nvindex;
	nIndex=0;

	fSinAngle=(float)sin(fAngle);
	fCosAngle=(float)cos(fAngle);
	fPos[3] = 1;
	for(ny=0;ny<para.nImageY;ny++)
	{
		fY=para.fCenterY+(ny-0.5f*para.nImageY+0.5f)*para.fSampleIntervalY;
		fPos[1] = fY;
		for(nx=0;nx<para.nImageX;nx++)
		{
			fX=para.fCenterX+(nx-0.5f*para.nImageX+0.5f)*para.fSampleIntervalX;
			fPos[0] = fX;
			
			fr=float(sqrt(fX*fX+fY*fY));
			
			if(para.bReconstructCylinder && fr > fFovRadius)
				continue;
			
			fU=fR+fX*fCosAngle+fY*fSinAngle;
			for(nz=0;nz<para.nImageZ;nz++)
			{
				fZ=para.fCenterZ+(nz-0.5f*para.nImageZ+0.5f)*para.fSampleIntervalZ;
				fPos[2] = fZ;

				ippmDotProduct_vav_32f(pfPmartix, sizeof(Ipp32f)*4,sizeof(Ipp32f), fPos, sizeof(Ipp32f), fResult, 4, 3);

				fnu=fResult[0]/fResult[2];
				
				nuindex=int(fnu);
				if(nuindex>=0 && nuindex<para.nDetectorColumnNumber-1)
				{
					fdeltau1=fnu-int(fnu);
					fdeltau2=1.0f-fdeltau1;
					fWeight=fCoefficient/(fU*fU);
					
					fnv=fResult[1]/fResult[2];
					nvindex=int(fnv);
					if(nvindex>=0 && nvindex<para.nDetectorRowNumber-1)
					{
						fdeltav1=fnv-int(fnv);
						fdeltav2=1.0f-fdeltav1;
						fui1=pfProjection[nvindex*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[nvindex*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						fui2=pfProjection[(nvindex+1)*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[(nvindex+1)*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						
						fuvi=fui1*fdeltav2+fui2*fdeltav1;
						pfImage[nz*para.nImageX*para.nImageY+ny*para.nImageX+nx]+=fWeight*fuvi;
					}
					else if(nvindex==para.nDetectorRowNumber-1)
					{
						fui1=pfProjection[nvindex*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[nvindex*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						pfImage[nz*para.nImageX*para.nImageY+ny*para.nImageX+nx]+=fWeight*fui1;
					}
				}
				
			}
		}
	}
}

void ComputeOneAngleEA(float *pPrj, float *pRcn, float fAngle)
{
	//////////////////////////////////////////////////////////////////////////
	//  filtering
	int nIndex=0;
	int nRow,nCol;
	float fDetectorX,fDetectorY,tmp;
	
	for(nRow=0;nRow<para.nDetectorRowNumber;nRow++)
	{
		nIndex=nRow*para.nDetectorColumnNumber;
		fDetectorY=(nRow+0.5f)*para.fDetectorTotalHeight/para.nDetectorRowNumber-0.5f*para.fDetectorTotalHeight;
		//////////////////////////////////////////////////////////////////////////
		//  Preweighting projection data
		tmp = para.fSourceDetectorDistance/sqrt(fDetectorY*fDetectorY+para.fSourceDetectorDistance*para.fSourceDetectorDistance);
		for(nCol=0;nCol<para.nDetectorColumnNumber;nCol++)
		{
			fDetectorX=(nCol+0.5f)*para.fDetectorTotalWidth/para.nDetectorColumnNumber-0.5f*para.fDetectorTotalWidth;
			pPrj[nIndex+nCol]=pPrj[nIndex+nCol]*cos(fDetectorX*PI/180)*tmp;
		}
		//////////////////////////////////////////////////////////////////////////
		//  convolution
		ippsConv_32f(pPrj+nIndex,para.nDetectorColumnNumber,pfConvKernal,para.nDetectorColumnNumber*2-1,pfBuffer);	
		memcpy(pPrj+nIndex,pfBuffer+para.nDetectorColumnNumber-1,para.nDetectorColumnNumber*sizeof(float));
	}
	
	//////////////////////////////////////////////////////////////////////////
	//  backprojection
	const float fdu=para.fDetectorTotalWidth/para.nDetectorColumnNumber;
	const float fdv=para.fDetectorTotalHeight/para.nDetectorRowNumber;
	const float fD=para.fSourceDetectorDistance;
	const float fR=para.fSourceAxisDistance;

	float fX,fY,fZ;	
	int nx,ny,nz;

	float fnu,fnv;
	float fSinAngle,fCosAngle;
	float fdeltaz;

	float fdeltau1,fdeltau2,fdeltav1,fdeltav2;
	float fui1,fui2,fuvi;
	float fr,fU,fWeight,fL2;

	const float fCoefficient=fR*fR*PI/para.nProjectionsPerRotation;
	const float fFovRadius=min(para.nImageX*para.fSampleIntervalX*0.5f,para.nImageY*para.fSampleIntervalY*0.5f);
	
	int nuindex,nvindex;
	nIndex=0;

	fSinAngle=(float)sin(fAngle);
	fCosAngle=(float)cos(fAngle);
	for(ny=0;ny<para.nImageY;ny++)
	{
		fY=para.fCenterY+(ny-0.5f*para.nImageY+0.5f)*para.fSampleIntervalY;
		for(nx=0;nx<para.nImageX;nx++)
		{
			fX=para.fCenterX+(nx-0.5f*para.nImageX+0.5f)*para.fSampleIntervalX;
			fr=float(sqrt((fX-para.fCenterX)*(fX-para.fCenterX)+(fY-para.fCenterY)*(fY-para.fCenterY)));

			if(para.bReconstructCylinder && fr > fFovRadius)
				continue;
			
			fU=fR+fX*fCosAngle+fY*fSinAngle;
			fnu=(atan((-fX*fSinAngle+fY*fCosAngle)/fU)*180.0/(fdu*PI))+para.fOffSet/fdu;
			nuindex=int(fnu+0.5f*(para.nDetectorColumnNumber-1));
			if(nuindex>0 && nuindex<para.nDetectorColumnNumber-1)
			{
				fdeltau1=fnu+(para.nDetectorColumnNumber-1)*0.5f-nuindex;
				fdeltau2=1.0f-fdeltau1;
				fL2=fU*fU+(-fX*fSinAngle+fY*fCosAngle)*(-fX*fSinAngle+fY*fCosAngle);
				fWeight=fCoefficient/fL2;
				for(nz=0;nz<para.nImageZ;nz++)
				{
					fZ=para.fCenterZ+(nz-0.5f*para.nImageZ+0.5f)*para.fSampleIntervalZ;
					fdeltaz=fZ-para.fStartSourceHeight;
					fnv=fD*fdeltaz/(fU*fdv);
					nvindex=int(fnv+0.5f*para.nDetectorRowNumber);
					if(nvindex>=0 && nvindex<para.nDetectorRowNumber-1)
					{
						fdeltav1=fnv+para.nDetectorRowNumber-int(fnv+para.nDetectorRowNumber);
						fdeltav2=1.0f-fdeltav1;
						fui1=pfProjection[nvindex*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[nvindex*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						fui2=pfProjection[(nvindex+1)*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[(nvindex+1)*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						
						fuvi=fui1*fdeltav2+fui2*fdeltav1;
						pfImage[nz*para.nImageX*para.nImageY+ny*para.nImageX+nx]+=fWeight*fuvi;
					}
					else if(nvindex==para.nDetectorRowNumber-1)
					{
						fui1=pfProjection[nvindex*para.nDetectorColumnNumber+nuindex]*fdeltau2+pfProjection[nvindex*para.nDetectorColumnNumber+nuindex+1]*fdeltau1;
						pfImage[nz*para.nImageX*para.nImageY+ny*para.nImageX+nx]+=fWeight*fui1;
					}
				}
			}

		}
	}
}
